home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / a2_8.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.5 KB  |  161 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 2.8 (Muller's Method).
  9. % Section    2.5,    Aitken's Process & Steffensen's & Muller's Methods, Page 97
  10. echo on; clc; format long; hold off; clear
  11. % This program implements Muller's method.
  12.  
  13. %    Define and store the function f(x) in the M-file  f.m
  14. % function y = f(x)
  15. % y = x.^3 - 3.*x + 2;
  16.  
  17. delete f.m
  18. diary f.m; disp('function y = f(x)');...
  19.            disp('y = x.^3 - 3.*x + 2;');...
  20. diary off;
  21.  
  22. % Remark. f.m and muller.m are used for Algorithm 2.8
  23. f(0); % Test for file f.m
  24. pause    % Press any key to see the graph y = f(x).
  25.  
  26. clc;
  27. %    Plot f(x) over the interval [a,b].
  28.  
  29. a = -2.5;
  30. b = 2.5;
  31. c = -5;
  32. d = 5;
  33. h = (b-a)/100;
  34. X = a:h:b;
  35. Y = f(X);
  36. axis([a b c d]);...
  37. plot(X,Y,'-g');...
  38. hold on;...
  39. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  40. xlabel('x');...
  41. ylabel('y');...
  42. title('Graph of y = f(x).');...
  43. grid;...
  44. axis;...
  45. hold off;...
  46. shg; pause    % Press any key to perform Muller's iteration.
  47.  
  48. clc;
  49. %    Place the starting values in  p0  p1  and  p2
  50.  
  51. %    Place the abscissa tolerance in  delta
  52.  
  53. %    Place the ordinate tolerance in  epsilon
  54.  
  55. %    Place the number of iterations in  max1
  56.  
  57. p0  = -2.6;
  58. p1  = -2.5;
  59. p2  = -2.4;
  60. delta = 1e-12;
  61. epsilon = 1e-12;
  62. max1 = 10;
  63.  
  64. [p2,y2,err,P] = muller('f',p0,p1,p2,delta,epsilon,max1);
  65.  
  66. pause    % Press any key for the list of iterations.
  67.  
  68. clc; clg;
  69. a = -2.65;
  70. b = -1.95;
  71. c = -9;
  72. d = 1;
  73. h = (b-a)/100;
  74. X = a:h:b;
  75. Y = f(X);
  76. max1 = length(P);
  77. n0 = min(7,max1);
  78. X0 = P(1:n0);
  79. Z0 = zeros(1,n0);
  80. axis([a b c d]);...
  81. plot(X,Y,'-g',X0,Z0,'or');...
  82. hold on;...
  83. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  84. xlabel('x');...
  85. ylabel('y');...
  86. title('Graphical analysis for Muller`s method.');...
  87. grid;...
  88. axis;...
  89. hold off;...
  90. shg; pause    % Press any key to continue.
  91.  
  92. J = 1:max1;
  93. Yp = f(P);
  94. points = [J;P;Yp];
  95. Mx1 = 'Iterations for Muller`s method.';
  96. Mx2 = ['     k                  p(k)               f(p(k))'];
  97. Mx3 = 'The solution is:';
  98. Mx4 = 'The error estimate for p is  ± ';
  99. clc,echo off, diary output,...
  100. disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
  101. disp('Iteration converged quadratically to the root.'),...
  102. disp(''),disp(Mx3),disp(''),disp('p = '),...
  103. disp(p2),disp('f(p) = '),disp(y2),...
  104. disp([Mx4,num2str(err)]),diary off,echo on
  105.  
  106. pause    % Press any key to perform Muller's iteration.
  107.  
  108. clc;
  109. %    Place the starting values in  p0  p1  and  p2
  110.  
  111. %    Place the abscissa tolerance in  delta
  112.  
  113. %    Place the ordinate tolerance in  epsilon
  114.  
  115. %    Place the number of iterations in  max1
  116.  
  117. p0  = 1.4;
  118. p1  = 1.3;
  119. p2  = 1.2;
  120. delta = 1e-12;
  121. epsilon = 1e-12;
  122. max1 = 10;
  123.  
  124. [p2,y2,err,P] = muller('f',p0,p1,p2,delta,epsilon,max1);
  125.  
  126. pause    % Press any key for the list of iterations.
  127.  
  128. clc; clg;
  129. a = 0.975;
  130. b = 1.425;
  131. c = -0.05;
  132. d = 0.55;
  133. h = (b-a)/100;
  134. X = a:h:b;
  135. Y = f(X);
  136. max1 = length(P);
  137. n0 = min(7,max1);
  138. X0 = P(1:n0);
  139. Z0 = zeros(1,n0);
  140. axis([a b c d]);...
  141. plot(X,Y,'-g',X0,Z0,'or');...
  142. hold on;...
  143. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  144. xlabel('x');...
  145. ylabel('y');...
  146. title('Graphical analysis for Muller`s method.');...
  147. grid;...
  148. axis;...
  149. hold off;...
  150. shg; pause    % Press any key to continue.
  151.  
  152. J = 1:max1;
  153. Yp = f(P);
  154. points = [J;P;Yp];
  155. clc,echo off, diary on,...
  156. disp(''), disp(Mx1),disp(''), disp(Mx2), disp(points'),...
  157. disp('Iteration converged quadratically to the root.'),...
  158. disp(''),disp(Mx3),disp(''),disp('p = '),...
  159. disp(p2),disp('f(p) = '),disp(y2),...
  160. disp([Mx4,num2str(err)]),diary off,echo on
  161.